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ABSTRACT 



An accurate method of parameter estimation for the mathemati- 
cal modeling of semiconductors using the Ebers-Moll equations is 
presented. Its usefulness is apparent in estimating parameters 
to be used in computer circuit-analysis programs that have been 
developed. The Ebers-Moll models were modified to better 
represent the actual characteristics. The least-square-error 
methods presented for estimating the parameters are easy to pro- 
gram to the digital computer and result in parameters that are 
quite accurate in describing the actual characteristics. This 
procedure yields solutions for parameters that are quite helpful 
to the engineer in solving electrical circuits involving p-n 
diodes and junction transistors. 



DUDLEY KNOX LIBRARY 
NAVAL POSTGRADUATE SCHOOL 
MONTEREY, CA 93943-5101 



TABLE OF CONTENTS 

I. INTRODUCTION 9 

II. EBERS-MOLL MODELS 13 

A. DIODE P-N JUNCTION MODEL 13 

B. JUNCTION TRANSISTOR MODEL 16 

III. LEAST-SQUARE-ERROR SOLUTION 24 

A. LEAS T-SQUARE -ERROR SOLUTION 24 

B. NORMALIZED LEAST-SQUARE-ERROR SOLUTION 26 

IV. DETERMINATION OF SEMICONDUCTOR PARAMETERS 29 

A. LINEAR EQUATION FOR DIODE PARAMETER ESTIMATION 30 

B. DETERMINATION OF PARAMETERS FOR JUNCTION 36 

TRANSISTORS 

C. TYPICAL VALUES FOR MODEL PARAMETERS 39 

D. COLLECTION OF DATA POINTS 42 

V. CONCLUSION 43 

APPENDIX Examples of Parameter Estimation 45 

COMPUTER PROGRAMS 56 

BIBLIOGRAPHY 61 

INITIAL DISTRIBUTION LIST 62 

FORM DD 1473 63 



3 



LIST OF TABLES 



I. Typical Values of Hp£^ 38 

II. Calculated Parameters and Standard Deviations 40 

III. A Comparison of Parameters Derived Using Two Methods — 41 



5 



LIST OF ILLUSTRATIONS 



Figure 

1. DC P-N Diode Circuit Model 

2. DC Junction Transistor Circuit Model 18 

3. Simplified Transistor Model 20 

4. Complete Diode Model 21 

5. Complete Transistor Model 23 

6. Plot of Measured Points for the 1N540 Diode SO 

7. Plot of Calculated Points Using the Least-Square 

Parameters for the 1NS40 Diode 51 

8 . Plot of Calculated Points Using the Normalized 

Least-Square Parameters for the 1N540 Diode 52 

9. Plot of the Measured Points for the 1N277 Diode 53 

10. Plot of Calculated Points Using the Least-Square 

Parameters for the 1N277 Diode 54 

11. Plot of the Calculated Points Using the 
Normalized Least-Square Parameters for the 1N277 

Diode 55 



7 



ACKNOWLEDGEMENT 



I wish to thank my thesis advisor. Professor Dr. S. R. Parker, 
for his assistance and advice . His guidance has helped me to 
investigate the material for this thesis. He has taken much 
of his time to discuss the problems that existed during this 
investigation . 



8 



I . INTRODUCTION 



In this age of computers, several excellent computer circuit- 
analysis programs have been developed which solve the equations 
describing electrical networks with remarkable accuracy. In order 
for the computed results to agree with actual circuit performance, 
the procedure of device modeling must furnish accurate models 
which are computationally efficient. As stated in a paper 
published by Design Automation, Inc. QlJ , "For all its advantages, 
the computer analysis of a circuit can be no better than the data 
that are fed into it. Poor models or improperly assigned values 
yield inaccurate results. In fact, once a circuit-analysis pro- 
gram is running on a computer and the designer has developed the 
skill to use it, the choice of semiconductor models and procedures 
for finding numerical values for the models’ parameters are the 
main stumbling blocks to applying the program. The parameters of 
a real switching transistor or real diode vary as functions of 
voltage, current and temperature. Modeling is most accurate when 
the device characteristics are measured under conditions of 
voltage, current and junction temperature similar to those under 
which the device is expected to operate." After a model has been 
chosen, it is essential that the various parameters must be 
computed easily and accurately. With nonlinear devices, such 
as diodes and transistors, the model equations must give a 
mathematical description of the curve so that for a given voltage 
the current is known or vice versa. For this analysis, the 
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Ebers-Moll [YJ models (modified slightly) for the diode and the 
transistor are usually employed. Using these models, methods 
are needed to determine the various parameters necessary for 
computation using any of the several computer circuit-analysis 
programs (TRAC, CIRCUS, NET I, NET II, and SCAN). These models 
cover both small-signal linear and large-signal nonlinear 
operation. The analysis in this thesis is limited to determining 
model parameters using steady-state characteristics. 

In May, 1965, engineers at North American Aviation, Inc. 
completed analysis on the least- square -error curve fit for the 
semiconductor diode model parameters for the SCAN and TRAC 
computer circuit analysis programs. The modified Ebers-Moll 
model for the diode was presented and the necessary assumptions 
were included. 

Design Automation, Inc., in June 1967, proposed a modeling 
procedure for transistors and diodes for computer-aided nonlinear 
circuit analysis [Y] . The models presented were essentially 
the Ebers-Moll models with additions. The analysis by Design 
Automations, Inc. was developed for the NET I Network Analysis 
Program. Several improvements were made to the existing models. 

In June 1968, the junction transistor model (modified Ebers- 
Moll model) was presented by E. Steele for the TRAC computer 
analysis program jjiJ . In March 1969 engineers at North American 
Rockwell Corporation presented, in an internal letter, a proposed 
method, using least-square estimation, of deriving the various 
parameters by analyzing the two junctions of the transistor 
separately . 
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In this thesis the problem of estimating the parameters of 
the Ebers-Moll models is considered. Questions investigated 
are: (1) determination of which regions of a semiconductor device 

are most useful for estimating each parameter of the model, 

(2) development of computer techniques of parameter estimation 
using least-square-error and normalized least square error 
solutions, and (3) estimation of the accuracy of the estimations. 
The ultimate question is whether the proposed models give an 
accurate representation of the actual performance of the semi- 
conductor devices. 

In Chapter TI the modified Ebers-Moll models for the diode 
and the transistor are presented. These models allow accurate 
calculations to be made of the semiconductor device parameters 
necessary for the computer circuit-analysis programs. 

In Chapter III the two methods of least-square-error 
solution are presented. One method is based upon the usual 
least-square-error solution and the other is based upon a normal- 
ized least-square-error solution as developed by Werther and 



Procedures for calculating the semiconductor device parameters 
are presented for both the diode and the junction transistor. A 
discussion of the useful regions of parameter determination and 
the typical values expected from this analysis is included. A 
method of accumulating measured data and correcting estimated 
values accordingly in also discussed. 

The parameters for typical diodes are calculated in the 
Appendix, using the two least-square-error methods. A comparison 
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of the results of these two methods is made and includes 
calculation of the error after the model parameters have been 
estimated. 

The computer programs used for the solution of the model 
parameters are given following the Appendix. Recursive equations 
for both least-square-error methods are used. The FORTRAN IV 
programming language is used for the computer programs. Most of 
the subroutines were developed for this analysis, but several 
were existing subroutines in the Naval Postgraduate School 
Computer Facility subroutine library. 
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II. EBERS-MOLL MODELS 



Ebers and Moll described the diode model and the junction 
transistor model, under DC conditions, using the small-signal 
parameters to relate to large -signal behavior [V]. These models 
described all regions of operation. Most of the other proposed 
models are mathematically equivalent and generally less accurate 



than the Ebers -Moll model L 5 J . Both the diode and the transistor 
models are considered below. The Ebers-Moll equations are modi- 
fied by considering a term called the emission constant M for the 
diode and Mg and Mg for the transistor. Also linear resistances 
are added to account for the deviations from the ideal models. 

The Ebers-Moll models relate the semiconductor device character- 
istics in terms of several parameters, and generally provide very 
little insight to the actual physical processes within the devices 
represented by these models . 



A. DIODE P-N JUNCTION MODEL 

For an ideal p-n junction diode, the current and voltage 
have the following relationship 

I = I s ' [exp@VAT) -T| (2.1) 

where I s is the reverse saturation current 
V is the voltage across the terminals 
q is the charge on an electron 
k is the Boltzman constant 

T is the absolute temperature (room temperature-298°K) 
kT/q equals 26mV at room temperature. 
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In practice very few p-n junctions can be described using the 
ideal diode relationship. In order to account for deviations 
between the theoretical and actual characteristics the ideal p-n 
junction diode needs to be modified by placing linear elements in 
parallel with or in series to the ideal diode. A large parallel 
resistance called the reverse diode leakage resistance is used to 
represent the slope of the reverse characteristics, and a small 
series resistance called the diode ohmic series resistance is used 
to represent part of the voltage drop across the p-n junction. 
Also, a modification to the Ebers-Moll model is the quantity 
M called the diode emission constant, with a value between 1 and 3, 
which changes the exponential dependence to more exactly describe 
the actual behavior, because of the non-ideal operation Qfj . 

The modified ideal diode equation becomes 

I = I s |exp(qV/MkT) - l] . (2.2) 

The DC diode circuit is given in Figure 1 using the above model. 
The various quantities included in Figure 1 are defined as 
Vj, Ij - the measured diode terminal voltage and current 
R g - the combined body, lead and contact resistance 
R-^ - the reverse surface-leakage resistance across the junction 
D - the modified ideal diode 

I - the ideal diode current given in equation 2.2. 

From the measured data of current and voltage, the parameters 

I , R , R and M are to be determined to describe a specific 
S s JL 

diode or p-n junction. 
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FIGURE 1 



DC P-N Diode Circuit Model 
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B. JUNCTION TRANSISTOR MODEL 

For the ideal junction transistor, Ebers and Moll proposed 
the following current-voltage relationship assuming a pnp 
transistor 



I P = I E0 



exp (qV^/kT) -1 



0( T 

I X CO 



1 - oCjjoCj 



exp (qV c AT) -1 



ic = 



= -^EO 






exp (qV^AT) -1 



L C0 






ixp (qV c AT) -1 



o( I o( I 

1 CO = N EO 



where 1^. = the emitter current 



( 2 . 3 ) 



( 2 . 4 ) 



( 2 . 5 ) 



I c = the collector current 
V£ = the emitter-to-base voltage 
Vq = the collector- to-base voltage 
^ = the common-base normal-mode DC current gain 
j = the common-base inverted-mode DC current gain 
Ieq = the saturation current of the emitter junction with 
zero collector current. 

Ic© = the saturation current of the collector junction with 
zero emitter current. 

From this the following is true QfJ 

^ES ~ I E0 // ^ 1_ °^N^l) ( 2 - 6 ) 



X CS I CO' / ^ 1_C ^N c (l) 



( 2 . 7 ) 



where I^ = the emitter-to-base p-n junction saturation current 



Iq 2 = the collector- to-base p-n junction saturation current. 
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In order to make the exponential dependence more exact, two 
quantities Mg and Mg are needed similar to the diode exponential 
emission constant. The quantity Mg is the emission constant for 
the emitter-base p-n junction and the quantity Mg is the emission 
constant for the collector-base p-n junction. The ideal transis- 
tor equations then become 



!g = I ES f xp (qVg/MgkT) -!]- O^Igg fxp (qV c /M c kT) -I] 



( 2 . 8 ) 



and 



Ir = 



= I cs [exp (qVg/MgkT) -I]- °<^- E s [e*P (qVg/MgkT) -I] . (2 . 9) 



Linear circuit elements are then added to modify the ideal transis- 
tor to an equivalent element to describe the actual relationship, 
similar to the diode formulation. The transistor equivalent 
circuit is presented in Figure 2 jjL^J . The elements of Figure 2 
are defined below: 

Ig = Ig + Ig the base current 

12 = the current generator across the collector-base junction 



= I 



cs fxp (qV c /M c kT) -I]' o( n I es [^xp(qV E /M E kT) -l] 

I-j^ = the current generator across the emitter-base junction 
= -*-ES f^*P (^Vg/MgkT) -lj- ^jlgg [^ X P (qVg/MgkT) -lj 
Rgg = the base spreading, bulk and contact resistance 
R ee = the emitter bulk and contact resistance 
R^ e = the collector bulk and contact resistance 
Rg = the emitter-base junction ohmic leakage resistance 
Rg = the collector- base junction ohmic leakage resistance 
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FIGURE 2 



DC Junction Transistor Circuit Model 
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The equivalent circuit for a pnp junction transistor may be used 
for a npn transistor by changing the polarities of Vp and Vp and 
the directions of Ij_ and I 2 . 

The specification sheets for transistors usually specify 

common-emitter DC current gains H FF and H FF rather than the 

tt n Jrt i 

common-base DC current gains and °( j. The relationship 

between these o( r s and Hpp’s are 

°^N = H FE n / ( H FE n +1 ) C 2 - 10 ) 

and 

°<I = H FE / (Hp E +1) . (2.11) 

1 1 

A simplified ideal model of the junction transistor is given 
in Figure 3, where Ip and 1 2 are replaced with diodes and 
dependent current sources . The resistances Rp and Rp are assumed 
to be very large and Rgg, Rpp, anc * k EE are assume( l to be small. 
For Figure 3 the elements are defined below 

I EF = I EsE Xp ( qV E /M E kT ^ "3 ( 2 - 12 ) 

and 

I CF = I C S E XP fe V c/Mc kT ) - 1 ] • ( 2 -13) 

From measured data, the parameters Rpp , Rpp, k BB’ Me, Mp, 

Ipg, Ips> q FE ’ anc ^ HpE . are to he determined to describe a 
specific transistor. 

The complete model including the AC components for the diode 
is given in Figure 4. The AC components added to the diode model 
are 

C-j- - the junction transition capacitance 
and 
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FIGURE 3 



Simplified Transistor Model 
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FIGURE 4 



Complete Diode Model 
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Cjj - the junction diffusion capacitance. 

In Figure 5 the AC components have been added to the 
transistor model so all operating conditions except the break- 
down regions for the junctions can be accounted in the model. 

The added components are 

Ct£ _ the emitter-base transition capacitance, 

- the collector-base transition capacitance, 

Cd£ - the emitter-base diffusion capacitance, 

and 

- the collector-base diffusion capacitance. 

The transition capacitances are functions of the junction voltages 
and the diffusion capacitances are functions of the emitted 
currents from the emitter and from the collector m. 
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FIGURE 5 



Complete Transistor Model 
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III. LEAST-SQUARE-ERROR SOLUTIONS 



The method chosen to compute the various parameters of the 
models presented is the least-square-error estimation. This 
solution is for a set of overdetermined simultaneous linear 
equations. There are two least-square-error methods presented. 
The first method is the usual least-square-error solution 
utilizing the concepts of the pseudo inverse as developed by 



the normalized least-square-error solution as developed by 



Werther and Parker LU A recursive equation is presented for 
each method to facilitate the use of additional data points in 
the calculation. These two methods and the estimations are 
compared with measurements to determine the usefulness of each 
me thod. 

The material for the least-square-error estimation and the 
normalized least-square-error estimation was obtained from 
Werther QT] . 

A. LEAST- SQUARE -ERROR SOLUTION 

The least-square-error solution is a very common method of 
solving a set of m simultaneous linear equations in n unknown, 
m being greater than n. This is the overdetermined case for 
linear equations which is common when measurements are made to 
find a small number of unknowns. The least-square- error solution 




Werther 




The second method is 
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solves the relationship"*" 

Ax = b (3.1) 

where A is the m by n matrix of coefficients, 
b is the m by 1 vector of constants, 
and 

x is the n by 1 vector of unknowns . 

A unique solution for x is obtained by using the concept 
of the Penrose pseudo inverse. This yields the exact solution 
if it exists, or the least-square-error approximation if the 
exact solution does not exist. The Penrose pseudo inverse is a 
n by n matrix and is identified as A + . The least-square-error 
solution is 




If the matrix A has rank n the least- square -error solution is 




T 

the matrix inverse of the product of A and A. In order for the 
rank of A to be n, the rows of A must be independent. This 
solution is easily obtained on the digital computer with available 
matrix operation subroutines. 

For the recursive relationship, the least-square -error 
solution for the first k points is 

% = p k A k T ^k C 3 - 4 ) 

where is equal to b^ which is the vector of k points, 

^The bar under a lower-case letter represents a column matrix 
or vector. 



25 



Aj< is the k by n coefficient matrix, 

T 

Pj^ is the inverse of the matrix A^, 



and 

'xk is the least-square -error solution. 

Now to consider the k+1 equation, a? is a row of new coefficients 
to the A matrix and is the new data for the vector b^+1- 



The re fore 



and 




(3.5) 



—k+1 



— k 



5k +1 



The k+1 least-square -error solution is 

= p ^ Tz 

-k+1 k+1 k+1 -k+1 

where Pj^+i is the inverse of the matrix ^Ak+l^k+l] 

The resultant recursive equations are 

A A Pj^a 

2ik+l “ *k + — f (zk+l-a T ^k) 

1+a Pj^a 



(3.6) 



(3.7) 



(3.8) 



and 



P k+1 " P k 



Pk-5. a T Pk 
l+a.^Pl<a 



(3.9) 



B. NORMALIZED LEAS T-SQUARE -ERROR SOLUTION 

The normalized least-square-error solution is a weighted 
least-square-error solution with weighting factors chosen such 
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that the solution lies as close as possible to all geometric loci 
described by equation 3.1. 

The matrix of weighting or normalizing factors is a diagonal 
m by m matrix called W. The factors are given by 



n 




3=1 



.,m 



(3.10) 



This weighting matrix is multiplied times the matrix A and 
vector b resulting in two new quantities 

A* = WA (3.11) 

and 



b* = Wb. (3.12) 

The normalized least-square -error solution is then 




(3.13) 



The recursive relationships for the normalized least-square- 
error solution are 



- £*. + p £ ** 

‘ k 1 l+a*p£ a*(z*k+i - £*x* k ) 



(3.14) 



and 



P* k+1 - P*k - 



p £ a*a* T pg 



l+a* T Pj a* 



(3.15) 



Using the fact that 



z* = (a T a) s z 
- k+1 ' “k+1 



(3.16) 
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and 



m 

a* = (a a) _a 

the recursive relationships become 



(3.17) 



A A p k — 

— *k+l = 2i*k + -m m" ' ■ (z k+1 - aTx*,) 

a^+a 1 ?^ a k 



and 



P* = p* 

k+1 k 



?k a a T P^ 



l+a T P* a 
— k — 



(3.18) 



(3.19) 



The recursive relationships for both the least-square-error 
solution and the normalized least-square-error solution are 
easily adapted to the digital computer and a sample computer pro- 
gram for each method is found after the Appendix. 
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IV. DETERMINATION OF SEMICONDUCTOR PARAMETERS 



Using the models presented in Chapter II, the model parameters 
are calculated utilizing the least-square-error method. Param- 
eter calculation is accomplished using the above equations and 
the equivalent circuits. In order to use the least-square -error 
methods, the nonlinear equations characteristic of these semi- 
conductor devices must be converted to a set of linear equations, 
which involves several assumptions. The solution is made on the 
nonlinear equations where these equations can be assumed to be 
linear. The error associated with these assumptions is so small 
it is considered negligible. 

These parameters enable the user to employ these models 
in existing computer circuit-analysis programs and obtain 
computed results that very accurately represent the actual per- 
formance. The data for parameter estimation may be obtained 
using the measured current and voltage values or values taken 
from device specification sheets. When specification sheets 
are used, the values of the voltage should be read to several 
millivolts in order to obtain accurate results. 

Since the data involves experimental error and is statistical 
in nature, as many data points as possible should be used in the 
estimation of the semiconductor parameters. Since the measure- 
ments are independent, in general, the rank of the matrix A in 
equation 3.1 is equal to n, the number of parameters. Since the 
over determined case must have the matrix A with a rank of at 
least n, the solution presented in Chapter III may be used. 



29 



A. LINEAR EQUATION FOR DIODE PARAMETER ESTIMATION 

The model for the semiconductor diode is presented in 
Figure 1. In order to determine the values of the parameters, 
the model equation has to be simplified and made linear to be 
compatible with the least-square-error methods. Using Figure 1, 
the following can be stated: 

I-j = I + V/R-l, (4.1) 

or 

I j = I s [exp(qV/MkT)-I] + V/R^ (4.2) 

The first restriction is that Vj (the measured terminal 
voltage) be positive (Vj^O). With this restriction, 
exp (qV/MkT) >> 1 because even with small positive measured 
voltages (for example 0.1 volt) and with q/MkT = 0.025, the term 
exp (qV/MkT) =54. 99 . This is much greater than 1. Since Rq is so 
very large compared to V, I V/R^, and equation 4.2 is simpli- 
fied to 

Ij = I s exp (qV/MkT) (4.3) 

which by rearranging becomes 

V = MkT/q In (I ./IJ . (4.4) 

Referring back to the model 

Vj = V + IjR s , (4.5) 

and substituting 4.4 into 4.5 results in 

Vj = MkT/q ln(Ij/l s ) + IjR s . (4.6) 
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This equation is simplified by letting Aq = MkT/q and 
Ajl = 1 /I S , obtaining 

Vj - Aq ln(A^Ij) + IjR s ( 4 . 7 ) 

V-j = Aq ln(A x ) + A 0 ln(Ij) + IjR s (4.8) 

Vj = A 3 + A 0 ln(Ij) + I jR s , (4.9) 

where A 3 = Aq ln(A^) . This approach is currently used by 
North American Rockwell Corp., among others. 

In the above equation, Vj and Ij are the measured forward- 
voltage and current data points for the diode. The other three 
quantities are constants to be determined using the least-square- 
error procedure. For the least-square-error method 

(4.10) 



(4.11) 

(4.12) 

where (x is the least-square-error solution. After A 3 , Ao and R s 
are determined, the parameters of primary interest are calculated 
using the following relationships: 



A = 



ln(Ii) 



and 



Vn 



b = 



V, 



m 



A 

x = 



a 3 

A 0 

Rc 
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M = A g qAT, 

A 1 = ex P (A 3 /A 0 ) , 



(4.13) 

(4.14) 



and 

I s = 1/A 1 . (4.15) 

The three parameters M, I s and R s are then calculated and 
can be used in the circuit-analysis computer programs. 

For determination of Ri, (the reverse diode leakage resist- 
ance due to the imperfections, surface effects, and thermally 
generated carriers in the depletion region) the same equation is 
used as the one for the determination of the other diode 
parameters. For this case, data points are used from the reverse 
characteristics of the diode ; thus Vj and I j are both negative . 

The basic equation is 4.2. By making an assumption that 
exp(qV/MkT) <<C 1 for V less than -0.2 volts, ^ equation 4.2 is 
reduced to 

q = - r s + (*+.lQ 

Solving for V, 

V “ (Ij + X S ) R 1 (4.17) 

and substituting this into equation 4.5 results in 

V j = (Ij + !s) R X + !j R s- C 1 *- 18 ) 

Rearranging, 

v j = J jC R i + R s ) + i s R i* ( 4 - 19 ) 

^With V = -0.2 volts and MkT/q = 0.026, exp (qV/MkT) = 4.5x10'^ 
and most measured data is less than - 0.2 volts. 
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Setting (R-^ + R s ) = Ag and IgR^ = Ag results in 

Vj = IjA 5 + Ag . (4.20) 

This is solved using the least-square-error procedures with 



A 





(4.21) 



(4.22) 

(4.23) 

In the above procedure Vj and Ij are the measured voltage 
and current in the reverse direction of the diode. From the 
results, R-^ is determined after an assumption is made. From the 
assumption that Rq>>R s , Rq is then equal to A 5 since A 5 = Rq + R s 
and R g is so small it can be neglected. A value of reverse 
saturation current, I s , may be found from this data but this 
value is not a good value for this current, especially for silicon 
diodes. For a 1N540 diode, the reverse saturation current is 
determined from the reverse characteristics along with Rq. The 
value obtained is then placed as a constant into the equation used 
to derive the forward parameters. The other two parameters are 
derived to see if this reverse saturation current resulted in 
parameters in the useable range. It is found that the value of 
R s is -1.46 ohms and the value of M is 3.67. These values are not 
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realistic numbers to be used in computer circuit-analysis pro- 
grams. These values are determined in the Appendix. The reverse 
characteristics of silicon diodes tend to be more difficult to 
model, since the ideal diode law breaks down completely as a 
description for the reverse silicon-diode behavior . In most 
p-n junctions the reverse current is not completely voltage- 
independent, due to the surface leakage. is the only reliable 

parameter derived from the reverse data. 

Several regions of the forward characteristics were investi- 
gated to discover where the best parameter estimations can be 
found. The standard deviation was computed to determine the 
closeness of fit to the measured points of the various regions. 
The closeness of fit is an indication of the error involved. The 
least standard deviation was found to be when the data was 
utilized along as much of the forward curve as possible . A 
comparison using the standard deviation from data received in 
different regions along the curve reveals this fact. Using the 
points corresponding to the large values of the measured voltage, 
or far out on the curve, the standard deviation is three times 
greater than the standard deviation found using the complete 
curve. The deviation found from the data of small values of 
voltage, or closest to the origin, is also much larger than the 
least standard deviation. These calculations are included in the 
Appendix using the 1N540 diode. 

The reverse saturation current determined from the forward 
characteristics was used to help find from the reverse 
characteristics by setting I s a constant in the derivation. This 
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method reveals that this value of I y has very little effect on 
the value of obtained. A simplified equation for determining 
RjL can be made, 

R 1 = V(l r l s ), (4.24) 

where Vj and Ij are measurements from the reverse characteristics. 

The parameters for two diodes, 1N540 and 1N277 are calculated 
in the Appendix. The standard deviation for each diode using the 
parameters determined is also given. The standard deviation was 
found using the measured voltage Vj and the calculated voltage VCj 
from the measured current I j . Equation 4.9 is solved placing the 
estimations of the parameters in the equation for the constants 
and using the various values of measured current. The value 
obtained is called VCj . The standard deviation is found by the 
following relationship 



standard deviation = 




(4.25) 



where N is the total number of data points used. This value 
of the standard deviation is used as an error indicator. The 
method with the largest standard deviation has the largest error 
using the same data points. For the two diodes used, it was 
found that the least-square-error method had a smaller standard 
deviation than the normalized least-square -error method, for 
both diodes. The difference is quite small and is not a good 
indicator of which method is best. From using these programs, 
it appears that the normalized least-square -error method has an 
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advantage in calculations when the measurements are at relatively 
great distances from the actual desired curve. Either method 
estimates parameters that give a curve that has less deviation 
than that expected due to variations from diode to diode due to 
manufacturing procedures and techniques. 

B. DETERMINATION OF PARAMETERS FOR JUNCTION TRANSISTORS 

The model for the junction transistor was presented in 
Figure 2. The parameters to be determined are Igg, I^g, Mg, ^C’ 
Rgg, ^cc > ^EE’ anc ^ R C - this analysis, the emitter-base and 

the collector-base junctions are considered to be two separate 
diodes for these calculations. For the emitter-base p-n junction 
parameters, measured data are made with the collector terminal 
open circuited and for the collector-base p-n junction parameters, 
measured data are made with the emitter terminal open circuited. 

Considering the emitter-base junction as a diode, calculations 
are made on the forward measured data to determine the quantities 
lES> Me, and RSE» the same as I s , M and R s were determined for the 
discrete diode above. The quantity Rgg is the combined emitter- 
base body, lead and contact resistance. The emitter reverse 
leakage resistance (Rg) is determined using the reverse character- 
istics, the same as for the diode. 

For the collector-base junction considered as a diode, 
calculations are made on the forward characteristics to determine 
the quantities Igg, ^C» and Rsc ’ the same as was computed for the 
diode. The quantity Rgg is the combined collector-base body. 
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lead and contact resistance. The collector reverse leakage 



resistance (R^) is determined using the reverse characteristics, 
the same as for the diode. 

To complete the model parameters for the junction transistor 
an assumption is made that the emitter bulk resistance is zero. 
The quantity is the emitter bulk resistance. Using this 

assumption, Rgg is calculated from 



R BB " ( H FE n + ^ R SE 



(4.26) 



The collector bulk resistance is determined to be 



R CC R SC “ R BB^^ H FEi + ^ (4.27) 

The quantities Hpp^ and Hpp^ are the common-emitter DC normal 

and inverted current gains respectively. 

The shift of H FF can be evaluated from the specification 
tfc n 

sheets available on the individual transistors. Specification 

sheets normally give curves of H FF versus I p . An average value 

irb n t 

is obtained by taking n values of Hpp from the curve at equally 

n 

spaced values of collector current over the range and then 
dividing by n. This average value of HpE n can then be used in 
further calculations of parameters. The parameter Hpp^ can be 
measured directly or assumed to be values as shown in Table I . 

The common-base 0('s can be calculated using equations 2.10 
and 2.11. 

This procedure of determining the model parameters for the 
junction transistor are currently used by North American 
Rockwell Corp . , among others . 
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TABLE I 

Typical Values of Hpp^ 

(Reprinted from the Motorola Switching Transistor Handbook.) 



Transistor Type Typical HpEi 

Planar, silicon epitaxial 0.2 

Alloy, uniform base 4 

Alloy, diffused base 1 

Mesa, diffused base 0.4 
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C. TYPICAL VALUES FOR MODEL PARAMETERS 

For the p-n junction diode the values of M, R and I are 

s s 

usually determined to be in the following ranges . The value of 

M is between 1.0 amd 3.0. the value of R s is between 10 “- 1 and 

10~1 kilohms, the value of I silicon is from 10“^-^ to 10”^mA and 

s 

for germanium is from 10 to 10 ^mA. The quantity R-^ has a value 
between one and ten megohms . 

The junction transistor parameters are usually contained in 
the following ranges: The value of Mg and Mg is between 1.0 

and 3.0, the value of Rgg and Rgg is between 10 and 

-2 - 4 i 

10 kilohms, and Rgg is usually between 10 and 10 kilohms. 

Rgg has been assumed to be zero in this analysis . The usual 

range of Igg is 10"-*- 1 * to 10"®mA for silicon and 10 - ^ to 10~^mA 

for germanium. The value of Igg is usually an order of magnitude 

larger than Igg . The quantities and Rgare usually larger 

than one megohm. 

These are only typical values, but if a parameter is 
substantially outside the usual range, the input data and 
calculations should be checked for error. These values were 
presented by Design Automation, Inc. frj. 

The parameter values for the two diodes calculated in the 
Appendix conform to the typical values for the model as is shown 
in Table II. 

In Table III a comparison of the values calculated using the 
least-square-error methods to the values calculated following 
the procedures given by Design Automation, Inc. There is 

very little difference in the parameter values derived by either 
method. The standard deviation for the least-square-error methods 
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M 

Is 

R s 

SD 



TABLE II 

Calculated Parameters and Standard Deviations 



1N540 



Usual Normalized 



1.75 

1.65xl0“ 10 

0.13 

0.006 



1.76 

1 . 86xl0 -10 
0.126 
0.0062 



1N277 



Usual 

1.07 

2.65xl0- 10 

82.83 

0.0025 



Normalized 

1.10 

3.72xl0" 10 
67 .9 
0.0027 
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TABLE III 

A Comparison of Parameters Derived 



Design Automation Least Square 



M 


1.75 


1, 


.75 


Is 


l.SxlO -10 


1. 


. 65x10”!° 


R s 


0.132 


0, 


.131 


SD 


0.0076 


0, 


.006 



Using Two Methods 

Normalized Least Square 
1.76 

1.86xl0- 10 

0.126 

0.0062 
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is a little less than that for the method of Design Automation, 
Inc. This indicates that the least-square-error methods are as 
well suited to these models as other procedures. 

D. COLLECTION OF DATA POINTS 

The measured values of current and voltage used in the 
A matrix and the Id vector for the least-square-error methods need 
to be evenly distributed from the minimum value of voltage (at 
least 0.1 volt) to a value of voltage that corresponds to maximum 
power. If these values are too close together or on only one 
portion of the curve, the derived parameters turn out to be 
negative values or values that do not correspond to the usual 
range. Most of this type of error is due to the inaccuracy of 
the measurements of the voltage and current. A typical derivation 
is presented in the Appendix. After the parameters have been 
initially estimated, the additional points are utilized in the 
recursive relationships . The additional points may be taken 
anywhere between the minimum and the maximum values of voltage . 

The least standard deviation is found utilizing as much of the 
forward characteristics as possible . 

For the reverse characteristics, fewer points are needed and 
they may be spread more than the forward characteristics . The 
reverse leakage resitance is a large value and if equation 4.24 
is used only one set of points is needed in the calculation. 

This calculation need not be too accurate since the reverse 
leakage resitance is of little importance . Values of current 
and voltage between -0.2 volts and the breakdown voltage are 
acceptable . 
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V. CONCLUSION 



A good model has been chosen since it has been shown that 
these models accurately represent electrical performance of the 
device over a very wide range and the models are computationally 
efficient in regards to the computer circuit-analysis programs 

f • ** 

presently in use. It is found that the error by this method is 
considerably less than the variation in the characteristics due 
to manufacturing techniques. The region chosen to measure the 
data for the parameters for the forward characteristics is 
represented as the points lying between 0.1 volts and the maximum 
power rating of the device. The measured data best represents 
the device when the points cover as much of the curve between the 
above mentioned points as possible. Considering the two methods 
of parameter approximation, both methods give results that are 
well within the accuracy expected. The normalized least-square- 
error solution does more to eliminate the influence of points 
located farthest from the desired curve. 

The standard deviation, which is a measure of the error in 
parameter estimation, reflects that these models give an accurate 
representation of the actual performance of the semiconductor 
devices. The standard deviation from the least-square-error 
methods is a few millivolts and most measurements could be this 
inaccurate. As seen on the forward characteristics presented in 
Figures 6, 7 and 8, the difference between the actual measurements 
and computed values is almost indistinguishable. 
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Parameters were not estimated for a transistor, but the 
method of determining these parameters is straight-forward and 
should not present any serious difficulty. Semiconductor devices 
with more p-n junctions than the transistor can be modeled 
similarly by dealing with each junction separately as a diode. 

A comparison of values of parameters from the least-square 
estimates with other methods of parameter estimation reveals that 
the least-square-error methods produce parameter values essential- 
ly the same as other procedures. The error for the least-square- 
error methods is slightly less than other procedures. 
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APPENDIX 



Examples of Parameter Estimation 

j 



This section is used to give two examples of the use of the 
procedures presented in this thesis. Two diodes, one the 1N540 

and the other the 1N277, were used in this analysis as the two 

{./ ' * 

examples. Measured data of current and voltage were taken in the 
forward and reverse portions of the two devices to be used to 

"5 • 

determine the device parameters. 

The following is the forward data collected for the two 

) . * ' l " 

devices with the voltage given in volts and the current in 
milliamperes . 

INS 40 1N277 



Voltage 


Current 


Voltage 


Current 


0.4 


0.002 


0.25 


0.002 


0.48 


0.006 


0.27 


0.004 


0.5 


0.0092 


0.28 


0.006 


0.52 


0.014 


0.285 


0.008 


0.54 


0.021 


0.29 


0.01 


0.56 


0.033 


0.3 


0.014 


0.58 


0.0555 


0.31 


0.02 


0.6 


0.086 


0.32 


0.025 


0.64 


0.215 


0.34 


0.049 


0.68 


0.483 


0.36 


0.085 


0.72 


1.15 


0.38 


0.14 


0.76 


3.18 


0.4 


0.222 


0.8 


9.15 


0.42 


0.37 


0.84 


15.5 






0.88 


37.8 






0.92 


76.0 






0.96 


155.0 






1.0 


270.0 







The reverse data for the devices are presented with both the 



current and the voltage as negative quantities and the current is 
in milliamperes and the voltage in volts . 
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IN 5 40 



1N277 



Voltage 


Current 


Voltage 


Current 


9.0 


0.01 


2.65 


0.004 


11.2 


0.012 


5.8 


0.008 


13.0 


0.014 


10.0 


0.0125 


15.5 


0.016 


11.5 


0.014 - 


19.5 


0.20 


20.0 


0.024 


22.5 


0.24 


30.0 


0.0355 


26.5 


0.0265 


36.0 


0.0425 



Using seven selected points, an A matrix is assembled to 
start the recursive operation with the forward characteristics. 
These seven selected points are spaced over as much of the curve 
as possible so as to take into consideration the influence of 
different parameters on various regions of the curve. Errors are 
sometimes encountered when two points are used that lie too close 
to each other. 

The seven points chosen for the A matrix and b vector are 
IN 5 40 1N277 



Voltage 


Current 


Voltage 


Current 


0.58 


0.0555 


0.25 


0.002 


0.6 


0.086 


0.27 


0.004 


0.68 


0.483 


0.285 


0.008 


0.72 


1.15 


0.3 


0.014 


0.88 


37.8 


0.34 


0.049 


0.96 


155.0 


0.38 


0.14 


1.0 


270.0 


0.42 


0.37 



where the voltage is in volts and the current is in milliamperes . 
The natural logarithms of the currents were determined and the 
information was then run on the computer programs to start the 
recursive procedures . The solution of the parameters using only 
the first seven points for the 1N540 are as follows for the 



least-square-error methods described above. 
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Least Square 



Normalized Least Square 



a 3 

A 0 

R s 

A 1 

Is 

M 



1.0249 

0.0453 

0.1268 

6.656xl0 9 

1.502x10-10 

1.743 



1.0233 

0.0451 

0.1324 

7 . 1373xl0 9 

1.4019x10-10 

1.734 



The solutions for the points for 1N277 are 

Least Square Normalized Least Square 



A3 


0.6097 


0.6162 


A0 


0.0275 


0.0281 


Rs 


78.2493 


71.3765 


Al 


4.1626xl0 9 


3.3997xl0 9 


Is 


2.4023x10-1° 


2.941x10-10 


M 


1.0587 


1.0798 



The rest of the measured points are used in the recursive 
equations and result in the final forward parameters as 

1N540 



Least Square Normalized Least Square 



A3 


1.0246 


1.0261 


A0 


0.0455 


0.0458 


Rs 


0.1308 


0.1257 


Al 


6.0536xl0 9 


5.3908xl0 9 


Is 


1.6519x10-1° 


1.855x10-1° 


M 


1.7495 


1.7611 



1N277 





Least Square 


Normalized Least Square 


A3 


0.6116 


0.6236 


A0 


0.0277 


0.0287 


Rs 


82.8309 


67.8968 


Al 


3.7768xl0 9 


2 . 6885xl0 9 


Is 


2.6477x10-1° 


3.7195x10-1° 


M 


1.0666 


1.1047 



The seven points from the reverse characteristics were only 
evaluated using the normalized least-square-error method since 
the only parameter attained from the reverse characteristics is 
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The value of is of little 



the reverse leakage resistance R^- 
importance and any value in the megohm range will be adequate.. 
The calculated values of R^ are « 

IN 5 40 1N277 

9.65xl0 5 8.69xl0 5 



Values of these parameters were calculated with the procedures 
of Design Automation, Inc . ,£ljusing the data as listed above. The 



following -values were obtained. 



1N540 



M 

Rs 

Is 



1.75 

0.132 



-10 



1.5x10' 

It is evident that these values are quite similar to those of the 
least-square-error methods. 

The standard deviation using equation 4.25 was calculated for 
the three sets of parameters calculated by the above methods. 

Each standard deviation calculated is listed below. 

Standard Deviations 
1N540 

Least Square 
0.006 



Design Automation 

. - !.: J'J. 0.0076 



Normalized Least Square 
0.0062 



1N277 

0.0025 0.0027 

In checking the various regions of the curve to be used the 
standard deviation was calculated for the parameters derived 
from data far out on the curve. The parameters calculated were 
M equals 1.872, I s equals 5.686xl0”10 and r s equals 0.101. The 
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standard deviation for these parameters is 0.018*+. This 
indicates that the best parameters are derived from data covering 
as much of the curve as possible . 

The standard deviation is a measure of the error between the 
calculated and the actual measurements. It is seen that the 
deviation is usually only a few milliamperes . 

Figures 6, 7, 8, 9, 10 and 11 are plots of the forward 
characteristics from the measured data and also of the calculated 
parameters for diodes. From these curves it is clear that the 
solution obtained by the least-square-error methods is very 
accurate . 
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FIGURE 6 



Plot of Measured Points for the 1N540 Diode 
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VOLTS 
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FIGURE 7 



Plot of Calculated Points Using the Least-Square 
Parameters for the 1N540 Mode 
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UNITS INCH. 
UNITS INCH. 



FIGURE 8 



Plot of Calculated Points Using the Normalized Least-Square 
Parameters for the 1N540 Diode 
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FIGURE 9 



Plot of the Measured Points for the 1N277 Diode 
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FIGURE 10 



Plot of Calculated Points Using the Least-Square 
Parameters for the 1N277 Diode 



54 



sc *wy 
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FIGURE 11 



Plot of the Calculated Points Using the Normalized 
Least-Square Parameters for the 1N277 Diode 
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COMPUTER PROGRAMS 



C THIS IS THE MAIN PROGRAM FOR THE NORMALIZED LEAST SQUARE 
C ERROR SOLUTION. 

C A IS THE M BY N INPUT MATRIX OF COEFFICIENTS. 

C B IS THE M BY 1 INPUT VECTOR OF CONSTANTS. 

C RS, SI, AND ETA ARE THE OUTPUT PARAMETERS WHERE RS 
C IS BODY RESISTANCE, SI IS THE REVERSE SATURATION CURRENT 
C AND ETA IS THE VALUE FOR THE EMISSION CONSTANT M.. 
IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION A (10,10) , B(10,10), C(10,10), D(10,10) 
DIMENSION W (10,10) , T(10), A2(10,10), E (10,10) 
DIMENSION A3 (10,10), F(10,10), G(10,10), H(10) 

READ (5 , 4) M,N 

READ (5 , 5) ( (A (I , J) , 1=1 , M) , J=1 , N 

READ (5 , 5) (B (1 ,1) ,I=1,M) 

4 FORMAT (2 110) 

5 FORMAT (7F10.0) 

- DATA T/10*0 . 0/ 

DO 118 1=1, M 
DO 118 J=1,N 
T (I) =T (I) +A (I , J) **2 
118 CONTINUE 

DO 122 1=1, M 
W (I , I) =1 . 0/DSQRT (T (I) ) 

122 CONTINUE 

DO 114 1=1, M 
B(I,1) =B (1 , 1) *W (I , I) 

DO 114 J=1,N 
A(I,J)=A(I,J)*W(I,I) 

114 CONTINUE 

DO 2*3 1=1, M 
DO 23 J=1,N 
A2 (J, I) =A(I , J) 

23 CONTINUE 

CALL MATMUL(A2 , N,M,A ,M, N, C ,M4,N4, L) 

CALL GAUSS3 (N,EPS ,C ,A3 ,KER, 10) 

IF (KER.EQ.2) GO TO 10 

CALL MATMUL(A3,M4,N4,A2,N,M,D,M5,N5,L2) 

CALL MATMUL(D,M5,N5,B,M,1,G,M8,N8,L3) 

PRINT 8, (G (1 , 1) , 1=1, N) 

8 FORMAT (10X, 3D14 . 6) 

DO 303 J=l,10 

READ (5 , 30) (H(I) ,I=1,N) 

READ (5, 32) Z 
30 FORMAT (3 F10.0) 

32 FORMAT (1F1 0.0) 

CALL WRLSF (A3,H,G,Z,N) 

PRINT 8, (G (1 , 1) , 1=1, N) 

A4=G (1, 1) 
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AO=G (2 , 1) 

RS=G (3 , 1) 

A1=DEXP (A 4/AO) 

PRINT 33, A4,A0,RS,A1 
33 FORMAT (10X, 4D14. 6) 

SI=1 . O/Al 
ETA=AO*38 . 46 
PRINT 8, SI, ETA 
303 CONTINUE 
GO TO 300 
10 WRITE (6, 12) 

12 FORMAT (10X,8HSINGULAR) 
300 STOP 
END 



C THIS IS THE MAIN PROGRAM FOR THE LEAST SQUARE ERROR 
C SOLUTION. 

C A IS THE M BY N INPUT MATRIX OF COEFFICIENTS. 

C B IS THE M BY 1 INPUT VECTOR OF CONSTANTS. 

C RS, SI, AND ETA ARE THE OUTPUT PARAMETERS WHERE RS 
C IS BODY RESISTANCE, SI IS THE REVERSE SATURATION CURRENT 
C AND ETA IS THE VALUE FOR THE EMISSION CONSTANT M. 
IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION A (10,10) , B(10,10), C(10,10), D(10,10) 
DIMENSION A2 (10,10) , A3 (10,10), G(10,10), H(10) 

READ (5,4) M,N 

READ (5,5) ((A(I,J) ,I=1,M) ,J=1,N) 

READ (5 , 5) (B(1 , 1) , 1=1, M) 

4 FORMAT (2 110) 

5 FORMAT (7F10.0) 

DO 23 1=1, M 
DO 23 J=1,N 
A2 (J , I) =A (I , J) 

23 CONTINUE 

CALL MATMUL (A2 , N,M, A ,M,N, C ,M4,N4, L) 

CALL GAUSS3 (N, EPS , C , A3 ,KER, 10) 

IF (KER.EQ.2) GO TO 10 

CALL MATMUL (A3 ,M4, N4, A2 ,N, M,D ,M5, N5 , L2) 

CALL MATMUL (D ,M5 , N5 , B , M, 1, G ,M8 , N8 , L3) 

PRINT 8, (G (1 , 1) , 1=1, N) 

8 FORMAT (10X, 3D14. 6 ) 

DO 303 J=l,10 

READ (5,30) (H(I) ,I=1,N) 

READ (5, 32) Z 
30 FORMAT (3 F 10.0) 

32 FORMAT (1F10.0) 

CALL TRLSF (A3 , H, G ,Z , N) 

PRINT 8, (G (1,1) , 1=1, N) 

A4=G (1,1) 

A0=G (2 , 1) 

RS=G (3 ,1) 

A1=DEXP (A 4/AO) 
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PRINT 33, A4,A0,RS,A1 
33 FORMAT (10X, 4D14. 6) 
SI=l. O/Al 
ETA=A0*38 . 46 
PRINT 8, SI, ETA 
303 CONTINUE 
GO TO 300 
10 WRITE (6,12) 

12 FORMAT (10X,8HSINGULAR) 
300 STOP 
END 



SUBROUTINE WRLSF (A3 , H, G ,Z ,M) 

IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION A3 (10,1) ,H(1) , Q(10) ,V2(10) , G(10,l) 
DIMENSION Q2 (10 , 10) , Q3 (10 , 10) , Gl (10) 

CALL MATVEC (A3 ,M,M,H,Q) 

CALL VECMUL(H,H,M,S) 

CALL VECMUL(H,Q,M,V) 

V1=S+V 

DO 122 1=1, M 
122 V2 (I) =Q (I) /VI 
DO 132 I-1,M 
Gl (I) =G (1,1) 

132 CONTINUE 

CALL VECMUL(H,Gl,M,Z2) 

Z3=Z-Z2 
DO 124 1=1, M 
12 4 V2 (I) =V2 (I) *Z3 
DO 126 1=1, M 
126 Gl (I) =Gl (I) +V2 (I) 

DO 134 1=1, M 
G (1 , 1) =Gl(I) 

134 CONTINUE 

CALL MULVEC (M,Q ,H,Q2) 

CALL MATMUL(Q2,M,M,A3,M,M,Q3,M,M,L9) 

Q4=l .ODO+V 
DO 128 1=1, M 
DO 128 J=1,M 

128 Q2 (I , J) = Q3 (I , J) /Q4 
DO 130 1=1, M 
DO 130 J=1,M 
A3 (I , J) =A3 (I,J) -Q2 (I , J) 

130 CONTINUE 
RETURN 
END 
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SUBROUTINE TRLSF (A3 ,H,G ,Z ,M) 

IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION A3 (10,1) ,H(1) , Q(10) , V2 (10) , G(10,l) 
DIMENSION Q2 (10,10) , <23(10,10) , Gl(10) 

CALL MATVEC (A3 ,M,M,H,Q) 

CALL VECMUL (H , Q , M , V) 

Vl=l . 0+V 
DO 122 1=1, M 
122 V2 (I) =Q (I) /VI 
DO 132 1=1, M 
G1 (I) =G(I,1) 

132 CONTINUE 

CALL VECMUL (H,G1,M,Z2) 

Z3=Z-Z2 
DO 124 1=1, M 
124 V2 (I) =V2 (I) *Z3 
DO 126 1=1, M 
126 G1 (I) =G1 (I) +V2 (I) 

DO 134 1=1, M 
G (1 , 1) =G1(I) 

134 CONTINUE 

CALL MULVEC (M,Q,H,Q2) 

CALL MATMUL (Q2 , M , M, A3, M , M , Q3 , M , M, L9) 

DO 128 1=1, M 
DO 128 J=1,M 
128 Q2(I,J)=Q3(I,J)/V1 
DO 130 1=1, M 
DO 130 J-1,M 
A3 (I , J) =A3 (I , J) -Q2 (I , J) 

130 CONTINUE 
RETURN 
END 



SUBROUTINE VECMUL (X, Y,M, SUM) 
IMPLICIT REAL*8 (A-H,0-Z) 
DIMENSION X(l) , Y(l) 

SUM = SUM-SUM 
DO S 1=1, M 
SUM=SUM+X (I) *Y (I) 

S CONTINUE 
RETURN 
END 



59 



SUBROUTINE MULVEC (M,X,Y,XY) 
IMPLICIT REAL*8 (A-H,0-Z) 
DIMENSION X(l) , Y(l) ,XY(10,1) 
DO 6 1=1, M 
DO 5 J=1,M 

5 XY (I , J) =X (I) *Y (J) 

6 CONTINUE 
RETURN 
END 



SUBROUTINE MATVEC (A , M, N , V , P) 
IMPLICIT REAL*8 (A-H,0-Z) 
DIMENSION A (10,1) ,V(1) ,P(1) 
DO 6 1=1, M 
P(I) =P(I)-P(I) 

DO 6 J=1,N 

P(I) =P (I) +A (I , J) *V ( J) 

6 CONTINUE 
RETURN 
END 



SUBROUTINE MATMUL(A,MA,NA,B ,MB ,NB, C ,MC ,NC ,L) 
IMPLICIT REAL8* (A-H,0-Z) 

DIMENSION A (10,1) ,B(10,1) ,C(10,1) 

L=0 

IF (NA.NE.MB) RETURN 

L=1 

MC=MA 

NC=NB 

DO 12 J=1,NC 
DO 12 1=1, MC 
C(I,J)=0.0 
DO 12 K=1,NA 

12 C(I,J)=C(I,J)+A(I,K)*B(K,J) 

RETURN 

END 



The other subroutine GAUSS3 is a subroutine in the Naval 
Postgraduate School Computer Facility subroutine library for 
calculating the inverse of a matrix if such an inverse exists. 

These programs are not the most computationally efficient but 
are adequate for this calculation. There are several improvements 
that can be made to these programs expecially in regard to 
storage required. 
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